Smectic and columnar ordering in length-polydisperse fluids of parallel hard cylinders 
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We apply a recently proposed density functional for mixtures of parallel hard cylinders, based 
on Rosenfeld's fundamental measure theory, to study the effect of length-polydispersity on the rel- 
ative stability between the smectic and columnar liquid crystal phases. To this purpose we derive 
from this functional an expression for the direct correlation function and use it to perform a bi- 
furcation analysis. We compare the results with those obtained with a second and a third virial 
approximation of this function. All three approximations lead to the same conclusion: there is a 
terminal polydispersity beyond which the smectic phase is less stable than the columnar phase. 
This result is in agreement with previous Monte Carlo simulations conducted on a freely rotating 
length-polydisperse hard spherocylinder fluid, although the theories always overestimate the termi- 
nal polydispersity because the nematic-columnar phase transition is first order and exhibits a wide 
coexistence gap. Both, the fundamental-measure functional and the third virial approximation, 
predict a metastable nematic-nematic demixing. Conversely, according to second virial approxima- 
tion this demixing might be stable at high values of the polydispersity, something that is observed 
neither in simulations nor in experiments. The results of the fundamental- measure functional are 
quantitatively superior to those obtained from the other two approximations. Thus this functional 
provides a promising route to map out the full phase diagram of this system. 



I. INTRODUCTION 

Simple fluids are made of atomic particles. These are 
identical, spherically symmetric particles which interact 
via a well defined interaction potential — of which the 
Lennard-Jones formula provides a very good approxi- 
mation. The classic states of matter are a consequence 
of this nature: interaction decays fast at long distance, 
hence we have a gas when the density is low; the potential 
has an attractive well at short distances, and this causes 
the liquid when the density is high enough; and finally, 
the interaction becomes strongly repulsive at very short 
distances and this makes the fluid freeze into a crystalline 
structure when it becomes very dense, due to entropic 
considerations (see e.g. jlj] for further details). 

In contrast to the 'simple fluid' paradigm provided by 
atomic fluids we have colloids. These are suspensions of 
big — around one micron — particles, which are actually 
aggregates of smaller particles, in a solvent which may 
also contain other elements (like ions, polymers, etc.). 
Because of this, particles are all different in shape, size, 
charge, etc., and the interactions are the result of adding 
up the separate contribution of each of the elements of 
the aggregate that we call a colloidal particle, as well as 
the entropic forces that the solvent and its constituents 
exert on them 0- Because of this, colloidal science has 
become the laboratory of liquid theory: almost anything 
in a colloid can be tuned so as to flt experimentally theo- 
retical models that would otherwise be considered highly 
academic (like hard spheres, to name the most famous 
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one). In particular, by varying the composition of the 
solvent one can modify entropic forces and gauge in this 
way the effective interactions between the colloidal par- 
ticles 0. 

One of the most interesting aspects of a colloid is its 
inherent polydispersity, i.e. the fact that colloidal parti- 
cles have different shape, charge, size, etc. The impact 
of this on the phase behaviour of the colloidal suspension 
is still a matter of active research. The study of polydis- 
persity is not new: Onsager himself paid attention to it 
in his famous 1949 article on the isotropic-nematic tran- 
sition of infinitely thin hard cylinders |^]. However, it is 
only in the last decade that the issue has regained the 
attention of liquid theorists, probably due to the funda- 
mental problems of formulating a statistical mechanics of 
polydisperse systems [5|. 

Most theoretical studies of polydisperse systems have 
focused on demixing and transitions between homoge- 
neous phases d 0, i 1 M, El, [11 m, Q El, [li. 

The reason is two-fold: on the one hand, specific tech- 
niques have been developed for those very common cases 
in which the free energy depends on the polydisper- 
sity distribution through a finite set of its moments 

[11,1131, or when polydispersity is small 
on the other hand, experimental data are available for 
colloidal liquid crystals and their transitions between the 
isotropic, nematic and non- uniform phases [23l[2^. How- 
ever, when it comes to finding theoretical approaches to 
spatial ordering transitions, results are more scarce due 
to the inherent difficulty of discerning how the ordering 
occurs in the continuum of species that form a polydis- 
perse mixture. In spite of this, several of these transitions 
have been tackled with different techniques. Interfaces 
and wetting have been successfully addressed with den- 
sity functional theory [H [H, [13, [M [H, [Ml , and so has 
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been freezing of polydisperse hard spheres, despite its 
higher difficulty [3l|, [33, [33| . 

Liquid-crystalHne spatially ordered phases (smectic 
and columnar) in polydisperse colloidal mixtures have 
received considerably less attention from the theoreti- 
cal point of view. Monodisperse fluids of hard rods are 
known to have a nematic-smectic transition 3J| . This is 
a continuous transition for parallel rods [35 1 which be- 
comes first-order when rods are allowed to freely rotate 
. In spite of some initial results that seemed to show a 
window of stability of the columnar phase [35| , it turned 
out that it was a finite-size effect, so that the smectic 
is more stable than the columnar for any aspect ratio of 
the rods [13 ■ It is also known that the addition of a sec- 
ond species of rods can destabilise the smectic phase in 
favour of a columnar phase [s^, [3^ . The same effect has 
been shown to occur in grand-canonical simulations of 
freely rotating, length-polydisperse, infinitely long rods 
(40| when polydispersity is larger than w 18%. There 
is also recent experimental confirmation of the enhance- 
ment of the stability of columnar ordering by polydisper- 
sity nil. 

The terminal polydispersity for the smectic phase had 
been predicted from a density-functional theory (4^ for 
a system of parallel hard cylinders. Despite the orienta- 
tional constraint, infinitely long rods are expected to be 
strongly aligned and thus behave very much as perfectly 
aligned rods — although not quite because the order of 
the nematic-smectic transition changes to first order for 
freely rotating rods, no matter their infinite length. The 
density functional used in [4^ was a very simple version 
of a weighted-density approximation (see e.g. [43] for a 
recent review), in which the weighting function is just 
proportional to the Mayer function. Simple as it might 
be, at the time there was no alternative density functional 
theory for carrying out this kind of analysis. But very 
recently a new functional for mixtures of parallel hard 
cylinders based on Rosenfeld's fundamental-measure the- 
ory 44 1 has been proposed [i^. The functional has been 
shown to provide an excellent estimate of the phase di- 
agram of the monodisperse system [37| . It is thus our 
aim in this paper to perform a bifurcation analysis of the 
smectic and columnar instabilities in this more accurate 
functional. 



plus an excess Fex[/o] = F[p] — i^id[p]- Here (3 = l/kT, 
with T the temperature and k the Boltzmann constant, 
and V(^) stands for the thermal volume of species I. The 
equilibrium density for a given chemical potential ^(l) is 
obtained from the Euler-Lagrange equation 



c(i)(r,0 = -(3 



■ exp 



m 

Sp{r,l)' 



{cW(r,0} 



(2) 
(3) 



If we specialise this equation for the uniform phase, of 
density profile ph{l), corresponding to the same chemical 
potential p{l)^ then 



PMO = ^exp{c«(0} 



(4) 



Here h(l) stands for the normalised probability density of 
particles of species I. From these two equation we obtain 



p(r, = exp {c« (r, - (0 } , 



(5) 



which will be the starting point of the bifurcation analy- 
sis. 



III. BIFURCATION ANALYSIS 

Let us assume that we have a length-polydisperse mix- 
ture of aligned hard cylinders in a nematic phase. A con- 
venient choice for I is I = L/ (L), where L is the length of 
the cylinders and (L) its average over the whole mixture. 
Let ph{l) be the density distribution of lengths in the ne- 
matic phase. Suppose that we reach a value of p at which 
the nematic fluid is no longer stable. Then the inhomoge- 
neous profile that emerges at the onset of the instability 
can be expressed as p{r, I) = ph(l) + e(r, I), where e(r, I) 
is a small perturbation. Using this expression in eq. ([5]) 
we obtain that, near the bifurcation point, 

e{r,l) = ph{l) I dl' I dr'c(r-r',/,r)e(r',r). (6) 



II. DENSITY FUNCTIONAL THEORY OF 
POLYDISPERSE MIXTURES 



where c(r - r', I, V) = -(35'^Fc^/5p(v, l)Sp{r', V) is the di- 
rect correlation function of the nematic phase. In Fourier 
space. 



Suppose that we have an inhomogeneous polydisperse 
mixture characterised by density profiles p{r,l), where I 
is a parameter (or set of parameters) which characterises 
the species (the length in our case). Then there is a 
Helmholtz free-energy functional F[p] which can be split 
into an ideal bit. 



PFi 



id 



[p] ^ j dlj dvp{v,l){\n[V{l)p{v,l)]-l], (1) 



e{ci,l) = ph(l) / drc(q,;,Oe(q,^'), 



(7) 



where as usual /(q) = / c?re*'i '"/(r). In order to pro- 
ceed we need to specify c(r — r',l,l') or, equivalently, 
i^ex- We will analyse three choices: (i) Fcx taken from 
the fundamental- measure density functional of ref. [4^ . 
(ii) a second virial approximation and (iii) a third virial 
approximation. 
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A. Fundamental-measure direct correlation 
function 

The expressions for the direct correlation function for 
this case appears as eq. (39) in ref. [i^ . For the case of a 
continuous polydisperse mixture, this is given, in Fourier 
space, by 

a, (3=0,1 

where q±^ and q\\ are, respectively, the lengths of the per- 
pendicular and parallel components of the wave vector in 
units of radius R and mean cylinder length (i) , 

c2;(")(g,0 = cos(qV2), (:>W(g,0 = ''"y^ (9) 



and V'ool'?) = 0, while 

+ y(i + 22/) 



q/2 



h 2yJo[q) 



^1(9) 



(10) 



^^2il^2y)Mq/-M 



2n / Ml) 



q 



(11) 



Here y = 77/(1 — 77), where rj — pttR^{L) is the total pack- 
ing fraction, and Jn{x) is the ri-th order Bessel functions 
of the first kind. 

The functional proposed in ref. [1^ is based on Tara- 
zona & Rosenfeld's functional for the fluid of hard disks 
dl]. We can use Rosenfeld's proposal for such a fluid 
instead [13] ■ Then the direct correlation function will 
still have the form ([5]), but the functions ipap{q) will be 
deflned as [4^] t/^qq ^ (q) = and 



4fiq)-i')'o'iq)^^y [Joiq)-Jiiq) 

+2(1 + 2y)Mq)^ + 277(1 + 2y) ^ '^'^'^^ 



/,(R) 



<^(9) = 22/2 
-2{l + 6y + 6y^ 



q \ q 

J^iq)-jUq) + 2{3 + Ay)Mq) 



Mq) 



(12) 

Jijq) 
q 

(13) 



If we now substitute ([8]) into (O, multiply the result 
by w(^)(q, I) and integrate over I, we obtain 

f7(q) =-E (E^7a(9||)V'a/3(<7x)) f^q), (M) 



where we have introduced the new functions 

C„(q) = y d/(:^(")(9||,0e(q,0, (15) 
TV^^gil) - I dlh{l)C:'^^\q\\M'-\q\\,l). (16) 



Eq. (|14p can be rewritten in matrix form as 

iJ(q)|(q) = [/ + iV(<7||) • vl/(q^)] ^(q) = 0, (17) 

where / is the 2x2 identity matrix, ^^(911) and ^'(g^) 
are the matrices with elements N^a{.q\\) and ipai3iq±) re- 
spectively and ^(q) is the vector with coordinates C/3(q)- 
Denoting 7i(q) = dct[i/(q)], the flrst non-trivial solu- 
tion of (fT7|) for which q 7^ follows from the couple of 
equations 



n{q) = 0, VH{q) ^ 0. 



(18) 



The flrst equation yields the value(s) of q for which 
^(q) 7^ 0, while the second one imposes that 7i(q) has 
a minimum at this value of q. These two equations de- 
termine the values of 77 and q at the bifurcation point. 
From dni), using dTO]), (HH) and HH), we obtain 

7i(q) = 1 + 2NQi{q\i)ijoi{q±) + Nn{q\\)i^ii{q±) 
+ [^01 (9|| ) -^00(911 )A^ii(9||)] ^^01(9^)- 

(19) 

Finally, the nematic-nematic dcmixing spinodal can be 
obtained as the value of rj at which 7i(0) — 0. 

Let us now consider three possible scenarios for a 
bifurcation: (i) nematic-nematic (N-N) demixing, (ii) 
nematic-smectic (N-Sm) bifurcation and (iii) nematic- 
columnar (N-C) bifurcation. 



1. N-N demixing 

Both versions, dTO]) and (HI]), and and ([13]), yield 
the same value for the functions 



V'oi(O) = 7/;io(0) =77(4 + 5^ + 277^), 
7/.n(0) = y\9 + Uy + 6y^). 

Also, from cqs. (fTB]) and ^ we have 

^oo(O) = A^oi(O) = 1, iVn(O) = l + A^, 



(20) 



(21) 



where A = y^Pj^^ is the standard deviation which 
characterises the degree of polydispersity of the system 
(remember that from our choice for I we have (l) = 1). 
We have introduced the short hand notation (/(/)) — 
J dh{l)f{l) for the mean value of a general function /(Z) 
with respect to the distribution fmiction h{l). Thus we 
flnd 



n{0) = 1 + 27/;oi(0) + 7/>n(0) + [7^11(0) - i'lii^)] A' 



(22) 
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The equation 7i(0) = leads to an analytic formula for 
the N-N demixing spinodal, namely 



'l + 2V7oi(0)+V.n(0) 



I 1 + 4:7] + 1]"^ 

7 — 2rj — Tj^ 



(23) 



Demixing appears for any rj > rj* where rj* is the solu- 
tion of (j23p for the maximum value of the polydispersity 
parameter A* (for the sake of reference, a Schultz distri- 
bution — see below — has A* = 1). 



2. N-Sm bifurcation 

A smectic instability is to be found by setting q± — 
and gjj = g > 0. Then introducing ([9]) in (|16p we obtain 



Noiiq) = iVio(g) = -(sin((zO), 

q 

Nn{q) = ^ [1 - (cos(90)] • 

q 



(24) 
(25) 
(26) 



Thus we obtain from ([191), & and p4)) - ((26|) . 
niq) = l + 2y(4 + 5y + 22/2)- 



2, (sin(gZ)) 



-I- 27/ 9 + 14y + 6y 5 

^ 2,. , . , o 2^2 [(sm(gO)^ + (cos(g/))^ - l] 

q^ 

(27) 

As a length-polydispersity model we make the stan- 
dard choice of a Schultz distribution function 

MO = ^^(^^^''exp[-(i.+ l);], z.>0, (28) 

whose mean is set to unity, i.e. (/) = / dllh{l) — 1. For 
this choice IS. — {u -\- 1)~^/^ < 1. It is easy to show that 
for this particular distribution function we obtain 

sin [A^^ arctan(gA^)] 



(sin(g/)) = 
(cos(gO> = 



cos [A~^ arctan((j'A^ 



[1 + g'A^] 



1/(2A2) 



(29) 



3. N-C bifurcation 

A columnar instability is to be found by setting = 
g > and q» = 0. Then using ((HI) in (HH) we find 



n{q) = 1 + 2^01(9) + V^ii(g) + [V>n(g) - 



Interestingly, ((30)) implies that the bifurcation to the 
columnar phase is independent of the particular func- 
tional form of h{l)^ it only depends on A, as it happens 
for the N-N demixing. 



B. Third and second virial approximations 

The third virial approximation of the direct correlation 
function for the system we are analysing has the expres- 
sion 



c(r,;,r) = /(r,;,r) 



l-f / dl"p{l")V{v,l,l\l") 



(31) 



where 



/(r, I, I') ^e{2R- r^) 9 {{I + l')/2 - \z\) , (32) 

with z expressed in units of (L), is the overlap function 
(minus the Mayer function) of two cylinders of the same 
radius R and reduced lengths I and I' . Q{x) is the Heav- 
iside step function (= if a; < and = 1 if x > 1). 
V{r,l,l' ,1"), the overlap volume between two cylinders 
of radius 2R and lengths I + I" and I' + I" , r being the 
vector joining their centres of mass, is given by 

V{rJ,l',l") = 8R^(L) (^arccosx - x\/l ~ a;2^ 

X e(i~x)x(^,^,z',n, (33) 

where x — r^/AR and 

x{z, I, r, I") = [{i + n/2 + 1" - \z\ - (|; - 1'\/2 - \z\) 

xQ{\l- l'\/2 - \z\)] 6 ({I + l')/2 + I" - \z\) . (34) 

The Fourier transform of — pc(r, ^, V) can be written in 
the same form ([5]), where now 



L q 

39 

^11(9) = -ri^U{2q\ 

TT 

U{q) = 16 / dxxM2qx) 



-C/(2g)77 



arccos x — x 



(35) 
(36) 



(37) 



and the functions uj^^\q,l) are given by eq. Thus, 
eq. (|19p becomes 



H(q) - 1 + 2^-01 (gi 



sin(g||0) 



■^ii{q±)' 



i^liiqi-) 



^l-(cos(gj|0)) 



((sin(g||0)^ + (cos(g||0)^ - 1) 



(30) 



(38) 
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We should point out that the second virial approximation 
can be obtained from (|36p just replacing U{2q) by zero 
(thus Vii(g) =0). 

The N-Sm {q± = 0) and N-C (qn = 0) bifurcations 
can be obtained from taking into account that 

U (0) in eq. (jST]) can be calculated analytically as U (0) = 
TT - 3V3/4. 

The uniform limit of ([55)1 yields 

7i(0) = 1 + 877 + 12cr?2 - 47^2 [4(1 + crjf - c] (39) 

where c = 4 — SVS/tt. Thus the spinodal of the N-N 
demixing is 



~ 277 y 4(1 + c?7)2 - c ' 



(40) 



The second virial approximation is obtained by setting 
c = in (|39p , which transforms the spinodal into 



Ab2 = 



477 



(41) 



IV. RESULTS 



The spinodals obtained from the second virial approx- 
imation of the direct correlation function are plotted in 
fig. [H As already mentioned, both the N-N demixing 
spinodal and the N-C spinodal are independent of the 
details of h[l), so they are valid for any polydisperse mix- 
ture. On the contrary, the N-Sm spinodal does depend 
on h{l). The curve of fig. [l]has been obtained using the 
Schultz distribution but in order to check what the 
effect of this choice is on this line we have also plotted 
the N-Sm spinodal for the distribution function 



r[(z.-i-i)/2]-+2 



T exp ■ 



r[{,y + 2)/2] 
r[(^ + i)/2] 



(42) 




which decays as a Gaussian and not as an exponential 
for long rods. For this choice 



(43) 



(and therefore A* = 0.755). The comparison of the N- 
Sm spinodal obtained with this distribution function with 
that obtained with the Schultz one reveals a weak de- 
pendence on h{l). For this reason we have stuck to the 
Schultz for the rest of the paper. 

The plot shows a crossover polydispersity. Ax = 0.394, 
below which the nematic bifurcates into a smectic and 
above which it does so into a columnar (from the h{l) 
given by (g^]) we obtain Ax = 0.428 instead). On 
the other hand, the N-N spinodal line reveals that N-N 



demixing can occur for very polydisperse mixtures (with 
A > 0.838). This is a defect of this approximation, as 
N-N demixing has never been observed in polydisperse 
systems of hard rods with a unimodal length distribu- 
tion. And not the only one, since an even more obvious 
drawback is the unphysical, high values of the packing 
fraction rj at which the spinodals appear. 

In striking contrast, the results provided by the 
fundamental-measure density functional proposed in 
ref. [45! (c.f. eqs. ^ and 1^, depicted 

in fig. mja), show a very different scenario. We also 
find a crossover polydispersity, at a slightly higher value 
Ax = 0.401. However, the N-N demixing is always 
metastable, as is consistent with simulations and exper- 
iments, and the values of the packing fraction at which 
the bifurcations occur are not far from the transition lines 
found in simulations. 

The same figure also shows the N-C spinodal resulting 
from the fundamental-measure density functional based 
on Rosenfeld's approximation for hard disks. It is most 
remarkable that for this functional no crossover is found. 
Thus this result either leads to the wrong conclusion 
that the smectic phase is more stable than the columnar 
phase for any polydisperse mixture, or it seems to sug- 
gest that the crossover polydispersity might be shifted to 
higher values. However, a definitive conclusion can only 
be achieved trough a coexistence calculation. 

Finally, fig. [IJb) shows the results obtained from the 
third virial approximation of the direct correlation func- 
tion (c.f. eqs. (US]) and ([551) 1. We can see a dramatic 
improvement with respect to the second virial approxi- 
mation in all details. N-N demixing becomes metastable 
and the values of the packing fraction at which the bi- 
furcations occur are much more reasonable. In fact, 
the scenario this approximation shows is rather close to 
the one obtained from the fundamental-measure density 
functional, the differences being only quantitative. 



V. DISCUSSION 

The phase behaviour of polydisperse mixtures of hard 
rods had received little theoretical attention mainly be- 
cause no good density functional theory was available for 
such a system, not even for the simplest model of aligned 
hard rods. Only very simple approximations, based on 
the Parsons-Lee rescaling, had been used. Despite the 
merit of these studies in finding a terminal polydisper- 
sity at which the N-Sm transition is preempted by a N-C 
one, this approximation is contingent on the accuracy of 
the second virial one — which we have seen not to be 
reliable for large polydispersity. 

Recently a functional based on Rosenfeld's fundamen- 
tal measure theory has been put forward for mixtures 
of parallel hard cylinders. In the present paper we 
have analysed its reliability in predicting the phase be- 
haviour of polydisperse mixtures. With this functional 
we also find a terminal polydispersity for the N-Sm tran- 
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FIG. 1: N-Sm and N-C spinodals obtained with the second 
virial approximation of the direct correlation function. The 
solid line always show the first phase (labelled in the figure) 
which bifurcates from the uniform nematic phase N. Length- 
polydispersity is taken as a Schultz distribution function, al- 
though only the N-Sm bifurcation line depends on this choice. 
The dot-dashed line shows the N-Sm spinodal obtained with 
the choice (|42|) for the polydispersity distribution. The filled 
circle and square show, respectively, the crossover polydisper- 
sity values Ax arising from this two distribution functions. 
Finally, the dotted line shows the location of the N-N demix- 
ing spinodal. 




sition and we confirm that N-N demixing can at best 
be metastable with respect to spatial ordering. We have 
also compared with the results obtained with second and 
third virial approximations. Although we find the former 
to have serious defects — like predicting N-N demixing 
at high polydispersity — the latter yields very reasonable 
results, close to those obtained with the fundamental- 
measure functional. 

Interestingly, a variant of the fundamental-measure 
functional constructed on Rosenfeld's proposal for the 
system of hard disk is not even able to predict the ter- 
minal polydispersity of the N-Sm transition. This calls 
for some caution in the use of Rosenfeld's functional to 
study the hard disk fluid. 

As for the validity of a bifurcation analysis, it obvi- 
ously provides the location of the phase transition if this 
is continuous, but it can be far from the coexistence line 
of the disordered phase for first order phase transitions. 
In the polydisperse system of hard rods, both the N-Sm 
and the N-C transitions are first order In the N- 

Sm transition this seems to be caused by the presence of 
particles aligned parallel to the smectic layers [i^. For 
this reason, in systems of perfectly aligned rods this tran- 
sition becomes continuous, so the N-Sm bifurcation line 
is the location of the transition predicted by the corre- 



FIG. 2: N-Sm and N-C spinodals obtained with the 
fundamental-measure density functional proposed in ref. [i^ 
(a) and with the third virial approximation (b). Lines and dot 
mean the same as in fig. [2] The dot-dashed line in (a) shows 
the N-C spinodal arising from the fundamental-measure den- 
sity functional arising from Rosenfeld's approximation for the 
hard disk fluid. Notice that it never meets the N-Sm bifurca- 
tion line. 



spending theories. The N-C is always found to exhibit a 
wide coexistence region both in simulations and in theory 
and therefore the terminal polydispersity found through 
a bifurcation analysis is but an upper bound of the true 
one. Polydispersity widens this coexistence region hence 
worsening the estimate provided by this bound. Locating 
this N-C coexistence is thus a necessary step to determine 
the N-Sm to N-C crossover. In the present state of the 
art this is a non-trivial task because the parallel hard 
cylinders functional contains a two-particle kernel which 
hinders the inclusion of polydispersity in inhomogeneous 
phases. How to circumvent this problem is a matter of 
current research. 
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